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Abstract — An unconditionally stable second order accurate implicit-implicit 
staggered procedure for the finite element solution of fully coupled thermoelas- 
ticity transient problems is proposed. The procedure is stabilized with a semi- 
algebraic augmentation technique. A comparative cost analysis reveals the su- 
periority of the proposed computational strategy to other conventional staggered 
procedures. Numerical examples of one and two-dimensional thermomechanical 
coupled problems demonstrate the accuracy of the proposed numerical solution 
algorithm. 


I. INTRODUCTION 

Transient response prediction of thermally loaded structures is of considerable 
importance in many aerospace engineering problems, and it has been the subject 
of intense research. Finite element formulations of the classical heat conduction 
problem without mechanical coupling have been presented by Wilson and Nickell 
[1]. Ritz type methods for the solution of linear dynamic problems in coupled 
thermoelasticity were given by Nickell and Sackman [2]. Oden [3] has formulated 
finite element models for the analysis of a class of nonlinear problems in dynamic 
coupled thermoelasticity, and Oden and Armstrong [4] have developed explicit 
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quadratic numerical schemes for the integration of nonlinear unpartitioned sys- 
tems of difference equations arising from the analysis of dynamic coupled ther- 
moviscoelastic problems. Recently, Ting and Chen [5] have introduced a unified 
numerical approach for the analysis of thermal stress waves. They have derived 
their algorithm from the concept of heat displacement and a variational formu- 
lation in Lagrangian form. They have proposed to integrate the resulting semi- 
discrete equations with conditionally stable explicit schemes. Liu and Zhang [6] 
have described an implicit-explicit procedure for the prediction of thermal stress 
waves in coupled thermoelasticity problems. They have adopted the explicit ra- 
tioned Runge-Kutta method [7, 8] for approximately solving the heat conduction 
equation and have claimed that their solution procedure is unconditionally sta- 
ble. However, their computational strategy requires the manipulation of a full 
matrix. In a sequel note, Liu and Chang [9] have slightly modified the original 
procedure of Liu and Zhang to involve a banded rather than full matrix, and have 
numerically verified the unconditional stability on one dimensional problems. 

However, several practical issues must be resolved before unconditionally sta- 
ble explicit rational Runge-Kutta schemes can become suitable for the analysis of 
real thermomechanical coupled problems. First, when unconditional stability is 
achieved for explicit time integration algorithms, typically consistency becomes 
conditional (see for example Hughes and Belytschko [10]). Second, most rational 
Runge-Kutta algorithms involve some divide operations by the difference be- 
tween intermediate solution quantities, which can significantly damage accuracy. 
Finally, these algorithms do not appear to accomodate staggered solution proce- 
dures for thermal/structure interaction problems, as they are not implemented 
in many existing production-level thermal computer programs. 

The semi-discrete equations governing soil-pore fluid interaction dynamic 
problems are similar to those governing thermoelastic coupled transient problems. 
In this sense, the work of Park [11], and very recently that of Zienkiewicz, Paul 
and Chan [12] can be extended to the response analysis of thermally loaded 
structures. 


In the present work, we present an unconditionally stable and robust implicit- 
implicit partitioned procedure for the solution of transient thermoelastic coupled 
problems. In Section II, we briefly review the basic equations for the linearized 
coupled thermoelasticity theory. A conventional implicit-implicit staggered solu- 
tion procedure is summarized in Section III. The thermal coupling term in the 
structural dynamics equation is treated as an applied force. However, while being 
very simple to implement, the resulting time integration algorithm suffers from 
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conditional stability. In Section IV, we introduce an augmented implicit-implicit 
staggered solution procedure for the partitioned problem. We establish the uncon- 
ditional stability and second order accuracy of the resulting numerical algorithm 
in Section V. In Section VI, we discuss the computer implementation aspects of 
the proposed computational strategy; we conduct a comparative cost analysis 
which demonstrates the superiority of the proposed solution procedure to other 
conventional staggered schemes. Finally in Section VII, we apply our partitioned 
algorithm to the solution of the one-dimensional Second Danilovskaya [13] and 
two-dimensional Youngdahl-Stemberg [14] problems. For both problems, the re- 
sults generated by the proposed stabilized procedure axe shown to be in excellent 
agreement with the analytical “exact” solutions. 


II. FINITE ELEMENT FORMULATION 


Let B denote the body of the structure to be analyzed, and dB = dB u U dB s U 
dBgUdBg the surface enclosing it. The basic equations for the linearized isotropic 
coupled thermoelasticity theory are: 


pii = diva- + b in B 

cO = —div( — kV9) — a(3A -{- 2fj.)8otr(e) + r in B 
a = 2/j.e + \(tre)I — a(3\ + 2p)(9 — 0 o )I 

e = i(Vu + Vu T ) 


and 

u = u on dB u 
an = s on dB s 
0 = 9 on dBe 
—IcVO = q on dB g 


( 1 ) 


where u, e, a , 6 , 9q, b, and r are the displacement, strain, stress, temperature, 
reference temperature chosen such that (9 — 9q)/9q << 1, body force, and heat 
supply fields, respectively, while p, A, c, a, p, k and n are the Lame’ moduli, the 
shear modulus, the specific heat, the coefficient of thermal expansion, the mass 
per unit volume, the thermal diffusivity, and the normal to the surface at a given 
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point, respectively. I is the identity tensor. The dot and T superscripts denote 
a time derivative and a transpose operation, and tr denotes the trace of a given 
tensor. 

If now we express the dependent variables u and 9 by suitable shape functions 
as: 


u = Nu and 9 = 

then a standard Galerkin procedure transforms 
coupled system of differential equations: 

Mii + Du + Ku — C9 
Q0 + H0 + 0 o C r u 

where M, D and K are the usual mass, damping and stiffness matrices, f is 
the prescribed structural loading vector, and Q, H, and r are respectively the 
capacity and conductivity matrices and the nodal source vector. If L denotes the 
differential operator corresponding to strain, the coupling matrix is expressed as 

C = / fi (LN) T (l, 1, 1, 0, 0, 0]Nd5. 

m. CONVENTIONAL IMPLICIT-IMPLICIT PROCEDURE 

In many applications, the coupling term C u that appears in the heat equation 
and which is induced by the effect of the strain rate is negligible. Therefore, 
one expects the second of equations (2) to remain parabolic and the temperature 
response to remain close to the uncoupled solution. Consequently, the dependent 
variable 9 is easier to predict than the displacement u, so that the most natural 
way of solving (2) would be: 

Mii n+1 + Du n+1 + Ku” +1 = f" +1 +C0 n+lP 

(3) 

Q0 n+1 +H0 n+1 = r n+1 - 0 o C T u n+1 


N9 

(1) in the following algebraic 
= f 

( 2 ) 


where # n+1 is the predicted temperature. Unfortunately, the above numerical 
procedure is only conditionally stable, even when each field is integrated with 
an unconditionally stable algorithm. Proofs of this result are given by Dubois- 
Pelerin [15] for various consistent predictors. Next, we introduce an augmentation 
technique that stabilizes the staggered solution of (2). 
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IV. AN AUGMENTED IMPLICIT-IMPLICIT PARTITIONED PRO- 
CEDURE 


Park, Felippa and DeRuntz [16] have introduced a differential augmentation con- 
cept that was successfully used in the stabilization of staggered solution pro- 
cedures for fluid-structure interaction problems. Basically, one of the coupled 
equations is injected into the other in order to “soften” the system, either by 
reducing the large eigenvalues of the uncoupled stiff equation, or by introducing 
some damping into it. Here, we adopt a different strategy. We perform a semi- 
algebraic augmentation — that is, we augment one of the two coupled equations 
while integrating both fields. 

First, the structural equation is integrated with the trapezoidal rule: 


u n+1 = ii n + =-(u n+1 -(- ii") 

= u n + ^(ii n + M _1 (r +1 - Du n+1 - Ku" +1 + C0 n+1 )] 

u" +1 = U n + y(u n+1 + u") 

At 2 

= u n + Atii" + [ii n + M -1 (F l+1 - Du n+1 - Ku n+1 + C0 n+1 )] 

4 ( 4 ) 


and the velocity vector is extracted as: 


(1 + -^M _1 D)ii n+1 = u n -(- -^-[u n -t- M -1 (f” +1 — Ku n+1 + C0 n+1 )] (5) 


Next, the heat equation is also integrated with the trapezoidal rule: 


d n+ 1 


A t , ■ n + l • n, 

= o n + — (d +e ) 

= e n + e n + Q‘'(r n+1 - H0 n+1 - ^ 0 C T u n+1 )] 


( 6 ) 


Finally, the system is augmented by recasting (5) in (6) to obtain: 
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( 7 ) 


9 n+ 1 = 9 n + ^-{9™ + Q -1 [r n+1 — H0 n+1 — + •^■M -I D) -1 

(u n + ^(ii n + M _1 (f ,+1 - Ku n+1 + C0" +1 )))]} 

Substituting (5) into the second of equations (4) and re-arranging (7) leads 


At 2 A i 2 

(1+ — B(At)M _1 K) u" +1 — B(At)M _1 C 9 n+1 = F n+1 

(--p^oAK) u n+1 + (1+ ^Q _1 H + ^0 o AC) 0 n+1 = R n+1 


(8) 


where 

A = Q _1 C t M -1 
B(A<) = (I + -^M _1 D)- ] 

F n+1 = u n + ^(I + B(At))u a + ^(B(A<)ii T, + M- , B(At)r +1 ) (9} 

V i 

R n +1 = Qn + + Q-I (r n+1 - Q 0 C T B(At)u n )] 

A/2 

- — [^ 0 Q _1 C T (B(AOu n + B(A<)M" 1 r+ 1 )] 


Now, a displacement predicted staggered procedure for the solution of (S) is: 

1. Predict the displacement field: 

u" +lP = u" (10) 

2. Solve for the temperature field: 

(I + ^Q-'H + ^qAC) 9 n+1 = R n+, + ^ 0 AKu" +lP (11) 
2 4 4 
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3. Correct the displacement field: 


A / 2 A / 2 

(1 + ^B(A0M _1 K) u n+1 = F n+1 + — B(A<)M _1 C 6 n+1 (12) 

4. Compute velocity, acceleration and flux fields: 

u n+1 = B(A<){u" + -^-[u n + M -1 (r* +1 — Ku n+1 + C0 n+1 )]} 

u n+1 = M _1 (r +1 +C0 n+1 -Du n+1 -Ku n+1 ) ( 13 ) 

0" +1 = Q-\r n + l -0 o C T u n+1 -H0 n+1 ) 

Remarks: 

p 

1. the predictor u n+1 is simply the previous step solution. It has been found 
(see Park [17]) that this is the most stable predictor when used in conjunction 
with the trapezoidal rule, while still maintaining a second-order accuracy. 

2. the injection of (5) into (6) is not arbitrary. It will be shown in Section VI 
that this is more economical than injecting (6) into (5). 

3. equations (13) define the computational path of the staggered procedure. 


V. STABILITY AND ACCURACY ANALYSES 


In this section, we establish that equations (10)-(13) result in a unconditionally 
stable second order accurate transient algorithm for the time integration of the 
coupled system (2). To avoid lengthy expressions, we consider the undamped 
(D = 0) and unforced (f = r = 0) case. Note however that even when D = 0, 
the quantity C9 still transmits a rate dependent damping effect to the structural 
equation. 


Stability. The stability of the proposed staggered procedure can be examined by 
seeking a nontrivial solution in the form: 


"u" +1 ' 
u n+1 

1+2 

U n ‘ 

u n 

•• n + 1 

n 

U 

— 

U 

1 

+ + 
c e 
<^> *<£> 
1 

1 — z 

e n 


( 14 ) 


i 


and determining under what condition the real part of z is positive. Substitution 
of (10) into (11) and (14) into (11)-(13) yields, after some algebraic manipulations: 


z 2 l + ^M -1 K -^pM _1 C 


u" 


'O' 

_ —(1 — z 2 )~0 o AK z 2 I + z% Q^H + ^oAC 


_e n . 


.0. 


(15) 


Therefore, the characteristic equation associated with (15) is: 


Mz 3 +VM^r 2 +(K+^ 0 CQ- 1 C T +— ^oCQ - 1 C t M _1 
2 4 


At 2 At 3 

K)— 2+VK 
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= 0 


(16) 


where 

V = cuc r 

U = Q -1 H(C r C) -1 


and | | denotes the matrix determinant. If the matrices M, K, Q and H are 
positive definite, and the coupling matrix C has full column rank, then U, V and 
each matrix coefficient of the determinant expression (16) is positive definite. If 
C is column rank deficient, U and V are positive semi-definite. In any case, all 
coefficients of the stability polynomial are non-negative. Consequently, the first 
part of the Routh-Hurwitz criterion [IS] for unconditional stability is satisfied. 
In order to check the second component of this criterion, we consider a 2-d.o.f. 
model problem for (2). The corresponding scalar form of (16) is: 


fl 3- 3 + Cl2 = ^ + a i * + oo = 0 


(17) 


where 


a 3 = 1 , a2 = 


Ath 
2 q ' 


a i = 


At 2 . 2 e 0 c 2 
-r\u + 




do = 


Afh 

Sq 


-c J 


Since At, h, q, u> 2 , 0 q, c 2 , and m > 0, then all the coefficients of the polynomial 
(17) in z are positive. Morevoer, the quantity 


Oi(Z2 — OqQ 3 


0 o hc 2 At 3 n 
Smq 2 1 + 



S 


is also positive, which demonstrates that the staggered solution procedure is un- 
conditionally stable for the 2-d.o.f. model problem. 

For the general multi-dimensional case, it turns out that the limiting case 
K = 0 which states that the structural system will grow quadratically in time, 
provides a sufficient test. For this case, (16) reduces to: 

| M 2 2 + jVMz + ■^•^ 0 CQ“ 1 C r | = 0 

Since M is positive definite and VM and CQ C are at least positive semi- 
definite, the procedure is unconditionally stable for the limiting case K = 0, as 
discussed in Bellman [19]. This argument has been extensively utilized in [12] 
during the analysis of several partitioned procedures. Therefore, we conclude 
that the procedure given by (10)-(13) is unconditionally stable. 

Remarks: 

1. the characteristic equation (16) reveals that the proposed procedure (10)- 
(13) is algorithmically identical to the one obtained by first differentiating 
the second of equations (2): 

Q0 + H0 + 0oC T ii = r 

then substituting ii from the first of equations (2) into the above equation: 
Q0 + H0 + 0 o C r M~ 1 C0 = = r - 0 o C T M~ 1 (f- Ku) 

However, differentiating the nodal source vector may be not practical, for 
example, if r is a discontinuous function of time. In our present derivation 
(11 )-( 13) we avoid this problem. 

2. the first-order thermal equation is algorithmically modified to behave as a 
damped second-order system. It should be emphasized that the described 
stabilization technique has not introduced any artificial damping. The only 
augmentation that is used is part of the governing equation of motion itself. 
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Accuracy. After differentiation, the third of equations (13) in the unforced case 
reads: 


9 


n+l 


= — 0 o Q -1 C T u n+1 - Q _1 H0 


-ll 


,n+l 


(18) 


Expanding the various terms in (8) around the time nAf and injecting (13) and 
(18) when needed leads to: 

Mii n + Ku n = C9 n + 0{At 2 ) 

Q9 n + H0 n = -6 0 C T u n + 0(At 2 ) (19) 


Comparing (2) and (19) demonstrates that the staggered procedure is second 
order accurate. The same result can be proved for the damped (D ^ 0) and 
forced (f ^ 0, r ^ 0) case. 


VI. COMPUTATIONAL ASPECTS 

In the remainder of this paper, we consider the case where the structure is un- 
damped (D = 0) and the mass and capacity matrices are lumped (M,Q are 
diagonal). The unconditionally stable staggered procedure (10)-(13) can be im- 
plemented as: 

1. Form: 


R n + !* = ^ r "+4Q[^ + ^(^_^ C V)]-^-^ C 7 '[ii n +M _1 (f +1 - K U n+lP )] 

( 20 ) 


2. Solve: 


(Q + ^H + ^^oC T M _1 C)<9 n+1 = R" +1 * 


( 21 ) 


3. Form: 


F n+1 * = M[u n + Af(u n + — ii n )l + — (f" +1 +C0 n+1 ) (22) 

4 4 ' 


4. Solve: 


10 


(23) 


5. Update: 


(M + 


A t 2 
4 


K)u" +1 


pn+l* 


u 


n+1 


U 


n+1 


0 


n + 1 


M" 1 (f l + 1 + C 0 


n+l 


— Ku n+1 ) 


u n + y(u n + ii n+1 ) 

Q-I( r n+1 _ ^ 0 c T u n+ i - H0" +1 ) 


(24) 


Equations (20) to (24) involve algebraic computations that are common to most 
implicit algorithms, when applied to the uncoupled problem. Only the quantity 
C~M -1 C deserves special attention. In particular, it is important to note that: 

• C M C is not a full matrix. It is a symmetric banded operator. Let n a , 
n*,, b 3 and bh denote the sizes and the semi-bandwidths of the structural and 
heat matrices, respectively. Typically, n, and b s are two to six times larger 
than rih and b^. The matrix product C M C is n/, by and heis a semi- 
bandwidth close to 2bh ■ Therefore, equation (21) entails the solution of an 
rift by rih symmetric banded system. On the other hand, if equation (6) had 
been injected into equation (5) — that is, if the temperature field had been 
eliminated from the structural equation — the resulting augmentation term 

1 T 1 . , 

would have been CQ _1 C' which is n, by n s and has a semi-bandwidth 
close to 2 bs. The latter would have entailed the solution of a symmetric 
system that is several times larger and denser than (21). For a rectilinear 
mesh composed of two-dimensional truss elements, the patterns of matrices 
C, C T , C t M - 1 C and CQ -1 C T are depicted in figure 1. 

• the additional cost incurred by the augmentation term is restricted to the 
factorization and subsequent solutions of equation (21). The precise value of 
this additional cost (with respect to the conventional procedure (3)) depends 
on the cleverness of the implementation. 
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I 


a 


FIG. 1 Patterns of the coupling matrices 
for a rectilinear mesh with 2D truss elements 
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In order to illustrate the computational costs of the proposed numerical 
procedure, we consider the problem of a clamped square plate where the edges 
are exposed to a sudden heating. The finite element mesh is composed of N by 
N 4-node regular elements. The stiffness and conductivity matrices K and Q are 
assumed to be stored in banded form so that operation counting is facilitated. 
In practice, these matrices are compacted in skyline data structures. We denote 
by d and p, respectively the number of structural degrees of freedom per node 
( d < 6) and the number of integration steps. 

The assumption of an N by N regular mesh with a number of fixed degrees 
of freedom at each node is unlikely in practice. However, it is the worst case as far 
as the computational effort required for the evaluation of the product CM C. 

For the above problem, the formation and factorization of equations (21) and 
(23) require respectively (2 + d)N 4 and d 3 N 4 / 2 multiplications. The resolution 
of equations (20)-(24) requires (7 cP + 6d + 3)N 3 multiplications for each time 
step. Therefore, the total computational effort needed for the transient coupled 
solution using the proposed stabilized procedure is: 

d 3 

E s ~ (— + d + 2)N 4 +p s (7d 2 +6d + 3)N 3 (25) 

For the same problem, the computational cost associated with a conventional 
second-order accurate conditionally stable procedure (3) is: 

E c ~ (^-^)iV 4 +p c {7d 2 + 6d + 3)N 3 (26) 


Clearly, unconditional stability is obtained at the cost of (d + 3/2)iV 4 addi- 
tional floating point operations. For linear problems, this computational effort is 
needed once. In the following, we show that this overhead is compensated by a 
much larger time step. 


The natural frequencies of the clamped square plate are given by: 


El 3 


U r 


= 7T 


m 2 + n 2 . 


12(1 — u 2 )p a 


(27) 


where E, v, l, and a are respectively Young modulus, Poisson’s coefficient, the 
plate thickness and its edge size [20]. Therefore, the lowest frequency is: 
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( 28 ) 


2n 2 I EP 

U7m,n - a 2 ]] 12(1 — u 2 )p 


and a good approximation of the highest element frequency is: 


EP 


J’t = 

m " o 2 V 12(1 - 1 A)p 


(29) 


An adequate time step for the stabilized procedure is given by u> m i n At 3 = 7r/10. 
For the conventional conditionally stable staggered procedure where both u and 
0 axe integrated with the trapezoidal rule, the stable time step is expressed as 
a multiple of the time step based on the Courant condition associated with the 
hyperbolic structural equation. Hence, A t c = m x 2/um\x, where m > 1. Using 
(28) and (29) we have: 


At 3 = 

“ 2 1/ 

f 12(l — u 2 )p 


A t c = 

20tt V 

ma 2 

EP 

/ 12(1 - u 2 )p 

(30) 

tt 2 N 2 

V EP 



so that 


p s = 40 

. 2i riV 2 (31) 

P = 

m 


axe the number of steps which would cover twice the largest period of the problem. 
The computational costs for both procedures become: 

E s ~ (^+d + 2)iV 4 

•V - (32) 

E c ~ —(7d 2 + 6d + 3)iV 5 

m 

which demonstrates the superiority of the proposed stabilized staggered procedure 
for N sufficiently large ( N > m/14). 
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VII. NUMERICAL EXAMPLES 


First, we consider the Second Danilovskaya problem [13]. An elastic half-space 
(x > 0) with the surface plane x = 0 assumed free of tractions for all time is 
exposed to a sudden high ambient temperature 9^. The continuum is assumed 
to be mechanically constrained and thermally insulated so that the displacement 
and temperature fields are given by: 


u x = u z (x,t) 
u y = 0 
u z = 0 
9 = 9(x,t ) 


(33) 


The boundary and initial conditions for this problem axe: 


<7**(0,<) = 0 


and 

u x (x,0) = 0 
it r (x,0) = 0 
0(x,O) = 9 0 


(34) 


where h is the boundary- layer conductance. The following dimensionless variables 
are introduced: 


ax 

K 

t 

K 

- ® IT 

Wo ( 35 ) 

9-9 0 
8 — 

_ _ a(A + 2^)u r 

k(39 0 
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where 


k_ 

pc 

a 2 _ A + 2/i 
• P 

0 = — - — 

3A + 2p 

The thermomechanical coupling parameter is defined by: 

s = th = Eh. 

pc{ A + 2/i) p 2 a 2 c 


( 36 ) 


(37) 


The exact solution for this problem can be obtained using the Laplace transform 
(see Nickell and Sackman [21]). The finite element solution is carried out using 
2-node linear elements. The ratio nh/ak is fixed to 0.5 and the thermomechanical 
coupling parameter S is set to 1. We report on the generated results for two time 
integration steps, At* 1 * = 7t/5 w miTl and At^ = Af^/2 = 7r/10 Wmin- These 
correspond to sampling the largest period of the mechanical problem into 10 
and 20 steps, respectively. Figure 2 depicts the dimensionless temperature 9 at 
x = 1.0 as a function of the dimensionless time t, for A t = At* 1 *. Figure 3 reports 
the dimensionless displacement u(f) at x = 1.0, for At = At* 2 *. As expected, the 
results for At = A are more accurate than those for At = A t^. However in 
both cases, the generated solutions are in good agreement with the exact ones. 
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FIG. 2 Dimensionless temperature at x = 1.0 
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FIG. 


3 Dimensionless displacement 


at x = 1.0 
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Next, we consider the case of an infinitely long elastic circular shaft of radius R , 
where the surface temperature undergoes a sudden uniform change over a finite 
band of length Z, and is steadily maintained thereafter (fig. 4). Youngdahl and 
Sternberg have presented in [14] an exact solution for the transient temperature 
and thermal stresses distributions in the shaft, when thermomechanical coupling 
is' neglected, in the form of definite integrals and infinite series. In cylindrical 
co-ordinates (r, <f > , z), the axisymmetric torsionless displacement and temperature 
fields axe given as: 


u r = u r (r,z,t) 
v.0 = 0 
Uz = u z (r,z,t) 
9 — 9(r,z,t ) 


(38) 


The boundary and initial conditions for this problem axe: 

a rr (R,z,t) — 0 
<r r z{R,z,t) = 0 

< 7 rT — » 0 as |z| — ► 00 
— * 0 as |z| -* 00 
o zz — * 0 as |z| — *■ 00 
c rz — »■ 0 as |z| — >■ 00 

o{R,z,t) = e 0 o |z| < ^ 

#(«,*. 0 = 0 w>| 


and 


u r (r,z,0) = 
u : (r, z,0) = 
Ur(r, 2 , 0 ) = 
U;(r,2,0) = 
0(r,z,O) = 


0 

0 

0 

0 

0 


(39) 


The following new dimensionless variables (there should be no confusion over the 
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present definition of these variables and their earlier use) axe introduced: 


P 

i 

6 

t 


R 

2z 

L 

e 

loo 

kt 


( 40 ) 


For ail computations, we set L = 2R and v = 1/4. The finite element 
solution is carried out using 4-node axisymmetric linear elements, and a time step 
At = 7r/10 u> mm . Figure 5 compares thepredicted temperatures at the center of 
the shaft ( p = 0) with the exact ones for 6 = 0, and reports on the effect of 
thermocoupling (6 = 0.5) on temperature distribution. Clearly, the stabilized 
procedure provides accurate solutions. The variations of the radial stress at 
£ = 0.1 for 6 = 0 and 6 = 0.5 are depicted in figure 6. All numerical results 
are reported at t = 0.2. It is interesting to note that when the thermocoupling 
effect is neglected the temperature field is overestimated, but the radial stress 
distribution is underestimated. 



FIG. 4 Problem geometry and finite element discretization 


20 




21 





FIG. 6 Dimensionless radial stress at f = 0.1 for t = 0.2 
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VIII. CONCLUSION 


An implicit-implicit staggered procedure for the solution of thermoelastic prob- 
lems is presented. It is stabilized with a cost-effective semi-algebraic augmen- 
tation scheme. The resulting transient algorithm is unconditionally stable and 
second-order accurate. 
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